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, We present an optimal strategy for the relative weighting of different data 

modes in inverse problems, and derive the maximum compatibility estimate 
(MCE) that corresponds to the maximum likelihood or maximum a posteriori 
estimates in the case of a single data mode. MCE is not explicitly dependent on 
the noise levels, scale factors or numbers of data points of the complementary 
i data modes, and can be determined without the mode weight parameters. As 

' a case study, we consider the problem of reconstructing the shape of a body in 

' from the boundary curves (profiles) and volumes (brightness values) of its 

\J^} , generalized projections. 
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1. Introduction 

In many inverse problems, various complementary data modes are avail- 
' able. For example, constructing the shape model of a body in R 3 is typically 

based on projectionlike data at various viewing geometries. In this paper, 
we consider the case where images I(u;,u;o) (generalized projections) ob- 
tained at viewing and illumination directions uj,luq E S 2 are available, but 
the reliable infomation in these images is only contained in the boundary 
curves dl between the dark background or a shadow and the illuminated 
portion of the target surface. This is a typical case in adaptive optics data in 
astrophysics, where the coverage of viewing geometries is also seldom wide 
enough to enable a full reconstruction of the model from images alone. 2 
Thus we include the possibility of augmenting the image dataset with a set 
of measured brightnesses (volumes of the generalized projections) L(uj,ujq) 
of the target at various observing geometries. 
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2. Case study: generalized projections 

We consider the inverse problem of determining the shape of a body BgR 3 
from some measured profiles of generalized projections dI(u>i,uJoi), i = 
1, . . . , n and their volumes L(oj oil w,), i = 1, . . . , m. 3 ' 4 

Our goal is to construct a total goodness-of-fit measure Xtot 

Xtot = xl + *dxl + *B.g(P), (i) 

where L denotes brightness data, d generalized profiles, and R regularizing 
functions g(P) (see Ref. 4 for discussion of these), where P G W is the 
vector of model parameters. Determining an optimal value for Xg (and Xr) 
is part of the inverse problem. 

The volumes of generalized projections are also called total or disk- 
integrated brightnesses: 3 

L(luq,lo) = / R(x; uj , u>)(ui, v{x)) da(x), (2) 

where A+ is the set of visible and illuminated points x G £>, 3 v(x) G S 2 
and da(x) are, respectively, the outward surface normal and surface patch 
of B, and i?(x;cj ,w) G R describes the intensity of scattered light at the 
point x on the surface. In its basic form, 

x \ =Y^[L(° hs \uJoi,LJi)-L(™ d \LJo i ,Ui)} 2 (3) 

i 

(assuming a constant noise level; see Ref. 3 and references therein for mod- 
ifications and variations of this), i-data on S 2 x S 2 uniquely determine a 
convex body and the solution is stable, 3 but L-data do not carry informa- 
tion on nonconvexities in most realistically available S 2 x S 2 geometries in 
practice. 

For many typical adaptive optics targets in our solar system, the gen- 
eralized profiles are starlike due to the proximity of uj and ujq and some 
regularity of the target shape at global scale. 2 Then we can write \% by 
considering, for each profile i, their observed and modelled maximal radii 
(from some point >f G M 2 within the profile) on the projection plane at 
direction angles ay (starting from a chosen coordinate direction): 

Xl = E^K-) - r^\a l0 )] 2 . (4) 

We now represent the body B as a polytope. Let two vertices a and b of 
a facet have projection points >c ai h^. The intersection point x of the radius 
line at a and the projection of the facet edge ab is readily determined. The 
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c (a) can now be determined by going through all eligible facet 
edges and their intersection points x a b(a): 

^ m ax d) («) = max { 1 1 x ab {a) - x 1| | a, b € V+ } , (5) 

where V+ is the set of vertices of the set of facets A+ approximating A+. The 
set A+ is determined by ray-tracing. 3 In general, facet edge circuits dA+ 
approximating dA+ (and corresponding forms of Xg) can be automatically 
derived for non-starlike shape models or profiles as well. 4 



3. Maximum compatibility estimate 

Let us choose as goodness-of-fit measures (from which probability distribu- 
tions can be constructed) the x 2 -functions of n data modes. Our task is to 
construct a joint Xtot with well-defined weighting for each data mode: 

n 

X 2 tot (P,D)=xi(P,D 1 )+J2^-iXi(P,D i ) D = {A,i = !,...,«} (6) 

i=2 

(to which regularization functions g(P) can be added), where D t denotes 
the data from the source i, and P G R p is the set of model parameter values. 
We assume the xf -space to be nondegenerate, i.e., 

arg min x f (P) ar g min x] ( p ) , * ¥= 3- 
In two dimensions, denote 

%W ■= {Xi|minxt ot ;A}, (7) 

yW ■= {x2l min Xto t ; A l- 

The curve 

5(A) := [log x(X), log y(\)} (8) 

resembles the well-known "L-curve" related to, e.g., Tikhonov regulariza- 
tion. ' However, here we make no assumptions on the shape of S. The 
curve S is a part of the boundary dlZ of the region KgK 2 formed by the 
mapping x ■ K p — > M. 2 from the parameter space P into xf-space: 

X = {P^ (logx?,log X 2)}, K = X {V) 

where the set V includes all the possible values of model parameters (as- 
suming that x is continuous and well-behaved such that a connected 1Z and 
dlZ exist). If the possible values of xf a rc n °t bounded, the remaining part 
dlZ \ S stretches droplet-like towards (oo, oo). The parameter A describes 
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the position on the interesting part S C dlZ, and it is up to us to define a 
criterion for choosing the optimal value of A. 

The logarithm ensures that the shape of 5(A) is invariant under unit or 
scale transforms in the x? as they merely translate S in the (logxf, logxl) - 
plane. It also provides a meaningful metric for the log xf -space: distances 
depict the relative difference in x 2 -sense, removing the problem of com- 
paring the absolute values of quite different types of xf • The endpoints of 
5(A) are at A — and A = oo, i.e., at the values of xf that result from 
using only one of the data modes in inversion. We can translate the origin 
of the (logXi, l°gX2)"Pl ane to a more natural position by choosing the new 
coordinate axes to pass through these endpoints. Denote 

x = loga;(A)|A = o = log minx? (9) 
Vo = log?/(A)|A-K>o = logminxa- 

Then the "ideal point" (xo,yo) is the new origin in the (\ogx, log y)-plane. 
A natural choice for an optimal location on S is the point closest to (xq, yo), 
i.e., the parameter values Pq G P such that 

P Q = argmin ([log X ?(P) - x } 2 + [logxi(P) - yo?) , (10) 

so we have, with A as argument, 

A = argmin ([loga;(A) - x } 2 + [logy(A) - yo} 2 ). (11) 

In this approach, neither the numbers of data points in each xf n or the 
noise levels as such affect the solution for the optimal Pq as their scaling 
effects cancel out in each quadratic term. Pq is thus a pure compatibility 
estimate describing the best model compromise explaining the datasets of 
different modes simultaneously. 

We call the point Po the maximum compatibility estimate (MCE), and Ao 
the maximum compatibility weight (MCW). This corresponds to the maxi- 
mum likelihood estimate in the case of one data mode, or to the maximum 
a posteriori estimate as well since we can include regularization functions 
here. If regularizing is used, the weights for the functions are either de- 
termined in a similar manner (see below), or they can be fixed and the 
regularization terms are absorbed in \\ (otherwise S C dlZ does not hold) . 

Another choice, frequently used in the L-curve approach, is to find the 
A at which S attains its maximum curvature, 1,6 but evaluating this point is 
less robust than finding Ao, and (11) is a more natural prescription, requir- 
ing no assumptions on the shape of S. We make two implicit assumptions 
here: 
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(1) The solutions Pqk corresponding to points on dTZ should be continuous 
(and one-to-one) in P-space along dTZ at least in the vicinity of the 
solution corresponding to A . If this is not true (in practice, if P\ = 
arg min x? t (-P) makes large jumps in P for various A around A ), one 
should be cautious about the uniqueness and stability of the chosen 
solution Prjj and restrict the regions of P included in the analysis. 

(2) The optimal point Ao on S should be feasible: if we have upper limits 
€i to acceptable %|, the feasible region T is the rectangle Hii^SXi < 
logej. If [logXi(-Po),logX2(-Po)] i J 7 and T C\TZ ^ 0, we choose the 
point on the portion S C TZ closest to the one corresponding to Ao (i.e., 
l°gX? = l°g £ i f° r one i)- H J 7 riTZ — 9, the data modes do not allow 
a compatible joint model, so either the model is incorrect for one or 
both data modes, or one or both e» have been estimated too low (e.g., 
systematic errors have not been taken into account). Note that model 
insufficiency should be taken into account in the estimation of . 

Note that, in the interpretation TZ = x(^); A, Xtot and dTZ are all in 
fact superfluous quantities, and we can locate the point estimate MCE Pq 
entirely without them with standard optimization procedures (and with no 
extra computational cost). However, it is useful (though computationally 
somewhat noisier) to approximate S via the minimization of Xtot 

with sam- 
ple values of A (see Fig. la), as in addition to obtaining the MCW Ao (and 
hence MCE as well) we can plot S to examine the mutual behaviour of the 
complementary data sources (including the position of the feasibility region 
T w.r.t. S). The solution for Ao is also needed for constructing distributions 
based on Xtot- Another possibility to examine TZ and dlZ is direct adaptive 
Monte Carlo sampling, but this is computationally slow. 

This approach straightforwardly generalizes to n x 2 -functions and n — 1 
parameters Aj describing the position on the n — 1-dimensional boundary 
surface dlZ of an n-dimensional domain TZ: the MCE is 

Po = argmin^[log^V i j , X % := minx? {?), (12) 
i=i x « 

and the MCW is, for A e E"" 1 , 

A = argmin^ [log * l ' tot } - 1 , xLot( A ) : = {xf\ minx t 2 ot ; a}. (13) 
Xio J 



i=l 



Another scale-invariant version of MCE can be constructed by plotting 



January 15, 2010 8:37 



WSPC - Proceedings Trim Size: 9in x 6in KaaslSAAC 



6 



Xi in units of Xi/x% ana ^ shifting the new origin to xf/xfo = 1 ; 



n 

Pq = arg min ^ |^ 




- 1 



A = arg min 



i=i 



n 



XiO 



1 . (14) 



This, however, is exactly the first-order approximation of (12) and (13) 
in S <C 1 when xl/xfo = 1 + $, giving virtually the same result as (12) and 
(13) as usually Xi(Po)/Xio — 1 <C 1 in the region around x? (Po), an d any 
larger ratios of Xi/xfo are n °t eligible for the optimal solution (see Fig. la). 

It is possible to use this approach for general regularizing functions g(P) 
as well (change x? — >■ g{P) for some i), but in such cases the shape of <S must 
be taken into account. If it is possible to have a solution g(P') = for a 
regularizing function g (or an almost vanishing g(P') such that log g(P') — > 
-co), one should, e.g., set a lower practical limit to g(P) by looking at the 
shape of 5, and choose the Ao within the restricted part of 5. 

4. Numerical implementation 

As examples of the optimal combining of brightness values and profile con- 
tours, we show some results for asteroid data. 

Fig. la depicts a typical evaluation of the curve S for 2 Pallas at various 
choices of A (or rather, this plot portrays the cross-section of the 2-surface 
dlZ in M 3 with smoothness regularization weight fixed at its final optimal 
value). The values for xj are normalized to be the rms deviations of model 
fits di — y/Xi/Ni, as in logarithmic scale this corresponds only to a shift 
of origin and a uniform linear change of plot scaling. The plotted points 
outline the curve 5(A) that is rather an oblique line than an L-shape, and 
the ideal point region, i.e., the point closest to the lower left-hand corner, 
can directly be found. The endpoints A = and A = oo stop at saturation 
regions rather than continue to large distances in the logx 2 -space. As can 
be seen from Fig. la, computational noise in the estimated points at A = 
and A = oo, corresponding to a small change of the position of the new 
origin w.r.t. S, does not affect the estimated location of the optimal point 
on S significantly. 

A sample observed vs. modelled profiles for 41 Daphne is shown in Fig. 
lb. The starlike surface model was described by the exponential Laplace 
(spherical harmonics) series for the surface radius r (for explicit positivity): 




(15) 



lm 
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Fig. 1. (a) <S curve plotted for 2 Pallas with various weights A (LC for brightness data, 
AO for adaptive optics profiles), (b) Sample observed (solid line) vs. modelled (dashed 
line) AO profile contour for 41 Daphne. Coordinates are in pixel units. 

truncated at I = 8, m = 6, with ci m as the shape parameters to be solved 
for. Other model parameters are the profile offset x for each image and the 
physical spin parameters describing the rotational transformations yielding 
the correct viewing and illumination directions (lj,lj ) on the body. 3-5 

5. Discussion 

The concept of the maximum compatibility estimate is directly applicable 
to any inverse problems with complementary data modes. The invariance 
properties of the MCE make it more generally usable than heuristic strate- 
gies for choosing the weights, especially when they use assumptions on the 
shape of dlZ or other case-specific characteristics. In our case study, the use 
of profiles is practical as it removes two sources of systematic errors inher- 
ent to using full images (brightness distributions X on the image plane): the 
errors in X from adaptive optics deconvolution and the model X errors due 
to the insufficently modcllablc light-scattering properties of the surface of 
the target body. 
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